UM-TH-99-02 
|hep-ph/ 9903468 
Spring Equinox, 1999 



WIMP Velocity Impact on Direct Dark Matter Searches 



Michal Brhlik a and Leszek Roszkowski 6 

a Randall Physics Lab 
University of Michigan 
Ann Arbor, MI 48109-1120, USA 

b Department of Physics 

Lancaster University 
Lancaster LAI 4YB, UK 



Abstract 

We examine the effect of some uncertainties in the input astrophysical param- 
eters on direct detection searches for WIMPs in the Galactic halo. We concentrate 
on the possible WIMP annual modulation signal recently reported by the DAMA 
Collaboration. We find that allowing for a reasonable uncertainty in a WIMP 
Maxwellian velocity distribution leads to significantly relaxed constraints on the 
WIMP mass as compared to the original DAMA analysis. 



1 Introduction 



It is generally assumed that spiral galaxies are surrounded by extended halos of dark 
matter (DM) for which a likely candidate is a hypothetical stable and neutral particle 
generically called a weakly interacting massive particle (WIMP). Experimental sensi- 
tivity has now improved to the level at which a signal from well-motivated WIMP 
candidates, like the neutralino of supersymmetry (SUSY), could be detected JT|. How- 
ever, inherent astrophysical and other uncertainties are likely to considerably weaken 
experimental limits and regions of a possible WIMP signal. 

One of the most promising ways to look for WIMPs in the halo of the Milky Way 
is the one based on annual modulation ||, [|. A WIMP signal rate depends on the 
effective velocity of an incident WIMP on a detector target which in turn depends on 
the Earth's velocity v e in the Galactic frame and on the WIMP velocity distribution in 
the halo. The Earth's motion in the Galactic frame is a superposition of its rotation 
around the Sun with that of the Sun around the Galactic center, 

(2Tt(t-t p )\ 

Ve [t) = V Q +V® COS 7 COS I I (1) 

where v Q ~ 232km/s and %, ~ 30km/s are the respective rotational velocities of the 
Sun and the Earth, 7 = 60° is the angle between the two planes of rotation, T = 1 yr 
and the peak occurs on t p = 153 rd day of the year (around the 2nd of June). Although 
the effect is expected to be small, around a few per cent, with enough sensitivity it is 
measurable. 

The annual modulation method has been adopted by the DAMA Collaboration 
in a Gran Sasso Laboratory-based experiment involving a detector consisting of nine 
9.70kg Nal(Tl) crystals. Based on the statistics of 14,962kgxday of data collected 
over a period from November '96 to July '97 (part of run II), the Collaboration has 
recently reported @] a statistically significant effect which could be caused by an annual 
modulation signal due to a WIMP with mass m x and WIMP-proton cross section a v 
given as 

m x = 59 GeVtH GeV £a p = 7.0^ x 1(T 6 pb (at 99.6% CL), (2) 

where £ = p x /po.3 stands for the local WIMP mass density p x normalized to po.3 = 
0.3GeV/cm 3 . (See also Figure 6 in Ref. [[| for a 2a signal region in the (m x ,£cr p ) 
plane.) According to DAMA, the new analysis is consistent with and confirms the 
Collaboration's earlier hint |J for the presence of the signal based on 4,549kgxday of 
data. 

In this Letter, we point out that the region selected by DAMA as corresponding 
to a possible WIMP signal strongly depends on the assumed values of astrophysical 
parameters which are poorly known. Given the fact that the Galactic halo has not been 
directly observed, nor any of its parameters have been actually measured, we advocate 
a more conservative approach. We demonstrate the effect by varying the position of 
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the peak of the WIMP velocity distribution which, following DAMA and a standard 
lore, we assume here to be Maxwellian. We find that significantly larger ranges of 
especially WIMP mass and, to a lesser extent, £cr p should be considered as consistent 
with the possible annual modulation effect reported by DAMA. This will have important 
implications for the allowed configurations of SUSY masses and couplings. 

In Section 2 we summarize the relevant elements of the DAMA analysis. In Sec- 
tion 3 we present the method which we use for calculating the expected event rate 
from Majorana WIMPs in annual modulation. Results and conclusions are presented 
in Section 4. 



2 DAMA Analysis 

The expected WIMP signal event rate (see Section 3 for details) is proportional to the 
elastic WIMP-nucleus cross section a(xN), the local WIMP number density p x /m x 
and the density of the target nuclei with mass tua- The local WIMP density enters 
multiplicatively and it is possible to account for its variation by normalizing it to some 
nominal value, e.g., po.3 m Ref. j§], and multiplying the scattering cross section by £. 
The event rate also depends on the velocity of a halo WIMP relative to the detector 
but this dependence is highly non-trivial and cannot be factorized out. 

Assuming the WIMP velocity distribution to be characterized by a parameter v , 
one can re-write Eq. (|T|) as 

V (t) =rj + Arj cos (^^j (3) 

where 77 = v e /v , r] = v Q /v , and Arj = v^cosj/vq. Since v@ 3> v® and it is expected 
that v ~ Vq, it is reasonable to assume that r] ^> Ar/. This was used by DAMA to 
approximate the differential signal event rate in the k-th energy bin by the first two 
terms in Taylor's expansion [|| 



S k [ V (t)} = S k [ Vo } + 



dSk 
drj 



A 77 cos 



'2ir(t - 1 



p) 



j 7/0 \ / 

= S , k + S m , k cos(^ Mt ~ tp) y (4) 

Using a maximum likelihood method @, [|], the (time-independent) background was 
separated from the constant component S 0jk and the time-dependent component S m>k . 
Because of the difference in the profiles of the expected energy distributions, the largest 
contribution to the signal was expected from the lowest energy bins above the thresh- 
old !• 

Based on the statistics of 14,962kgxday of data, values of S 0>k and S m ^ k were ob- 
tained between 2keV (the software threshold) and 20keV as shown in the Table @]. 
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Energy 
(keV) 


Sq^j~ ± &o^k 
(cpd/kg/keV) 


(cpd/kg/keV) 


2 - 3 


0.54 ±0.15 


0.018 ±0.009 


3 - 4 


0.23 ±0.08 


0.012 ±0.004 


4 - 5 


0.09 ±0.04 


0.006 ±0.002 


5 - 6 


0.04 ±0.02 


0.003 ±0.001 



Above 6 keV the detected time-dependent component is statistically insignificant. 

As mentioned above, it is possible to write So^ = £cr p S ok and S m ^ = i^ p S mk ||. 
Both S'q k and fe depend on WIMP mass but also on its velocity distribution, in 
addition to several other detector and nuclear factors. Expressions for S' k and S m k 
used by DAMA can be found in Ref. ||. The ranges (0) of m x and £a p in agreement 
with the hypothesis of an annual modulation signal due to a Majorana WIMP with 
dominant spin independent interaction with Nal target were obtained by a maximum 
likelihood method [f| . The period was assumed to be one year and the position of the 
peak was fixed at the expected value. 

We note that in DAMA's analysis a spherical halo model and a Maxwellian velocity 
distribution with a fixed value of vq = J 2/ 3v ims = 220km/s was assumed [Q. We will 
argue that the ranges of m x and £cr p depend sensitively on the assumed value of t>o- 
First we will present the procedure which we use. 

3 Majorana WIMP Detection Rate in Nal Detec- 
tors 

Elastic interactions of relic WIMPs with nuclei in the detector result in a nuclear recoil 
energy deposit in the detector volume. Their size depends on the cross section of 
the WIMP scattering off constituent quarks and gluons. For non-relativistic Majorana 
particles, these can be divided into two separate types JFj . The coherent part described 
by an effective scalar coupling between the WIMP and the nucleus is proportional to the 
number of nucleons in the nucleus. It receives a tree-level contribution from scattering 
off quarks, xi ~~ > X<2S as described by a Lagrangian C ~ (xx) (<?<?)• The incoherent 
component of the WIMP-nucleus cross section results from an axial current interaction 
of a WIMP with constituent quarks, given by C ~ (xi^lnX) (97^759); an d couples the 
spin of the WIMP to the total spin of the nucleus. 

In the case of a supersymmetric neutralino, there are several diagrams contributing 
to the scalar part, mainly from Higgs exchange and squark exchange with, in general, 
non-degenerate left and right squark soft masses. Another contribution to coherent 
interactions comes from one-loop neutralino-gluon scattering where the exchanged Higgs 
couples to a heavy quark loop on a gluon line [Q . The axial interaction is due to the Z 
and squark exchange. 

The differential cross section for a WIMP scattering off a nucleus X% with mass rriA 
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is therefore expressed as 



da da scalar da axial 



d\q\ 2 d\q} 2 d\q\< 

where the transferred momentum q = mAm-x v depends on the velocity v of the incident 

WIMP. The effective WIMP-nucleon cross sections <j scalar and <j axial are computed by 
evaluating nucleonic matrix elements of corresponding WIMP-quark and WIMP-gluon 
interaction operators. 

In the scalar part contributions from individual nucleons in the nucleus add coher- 
ently and the finite size effects are accounted for by including the scalar nuclear form 
factor F(q). (The effective interaction in general also includes tensor components but 
the relevant nucleonic matrix elements can be expanded in the low momentum-transfer 
limit in terms of the nucleon four-momentum and the quark (gluon) parton distribu- 
tion function. As a result, the non-relativistic WIMP-nucleon Lagrangian contains only 
scalar interaction terms.) The differential cross section for the scalar part then takes 
the form |l 

J scalar i 

— = - [Zh+[ A-Z)UV {q ), (6) 

where f p and f n are the effective neutralino couplings to protons and neutrons, respec- 
tively. Explicit expressions for the case of the supersymmetric neutralino can be found 
in the Appendix of Ref . |§ . 

The effective axial WIMP coupling to the nucleus depends on the spin content of 
the nucleon Ag p n and the overall expectation value of the nucleon group spin in the 
nucleus < S p>n >. For a nucleus with a total angular momentum J we have 

J^axial o 

■#=^ + D%), (7) 
with A = j[a P {S p ) + a n {S n )}. The axial couplings 

1 ^ ' ' (p) n - _L A A» (9\ 



V * u d s V ^ uA* 



are determined by the experimental values of the spin constants Au^ = Ad^ = 0.78, 
Ad (p) = Au in) = -0.5 and As (p) = As (n) = -0.16. The effective couplings d q depend 
on the WIMP properties and for the neutralino they can be found in the Appendix of 
Ref. !• 

In translating <j(xq_) and cr(xg) into the WIMP-nucleon cross section in Eq. (f|) 
several uncertainties arise. The nucleonic matrix element coefficients for the scalar 
interaction are not precisely known. Also, the spin content of the nucleon and the 
expectation values of the proton (neutron) group spin in a particular nucleus are fraught 
with significant uncertainty and nuclear model dependence. These ambiguities have to 
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be considered in numerical calculations. Finally, in order to obtain a(xN), models of 
nuclear wave functions must be used. The scalar nuclear form factor reflects the mass 
density density distribution in the nucleus. Following DAMA analysis we take the form 
factor for Na to be equal to one while for I we use the Saxon- Woods form factor [|Kj 



F(q) = ^Jl^A e ^)\ (9) 
qRi 

where R\ = V R 2 — 5s 2 , R = x 1.2 fm, j 1 is a spherical Bessel function and s = 1 fm. 
The spin form factor for I is assumed to be 

S(q) = [0.6845e- 359 - lg2 - 3.427g 2 + 0.30A2}e~°- 2q2{r "- r2 ^ ) (10) 

where q is in GeV and the nuclear radius is r = l^A^Hm. The finite size effects for Na 
can be neglected. The nucleon group spin expectation values for I have been estimated 
within the quenched interacting boson model giving < S„ >= 0.128 and < S n >= 



TT| , and for Na we use the odd group model result < S p >= 0.136 and < S n >= [12 

In order to calculate the WIMP detection rate for a given material, it is necessary to 
convolute the cross section Eq. (|5|) with the local WIMP flux which, in the rest frame of 
the detector, will be time-dependent. The differential event rate can be expressed as [[p] 



T(Ef>)< [Zf p + (A — Z) f n ] F (q) 
+8A 2 J(J+l)S(q)\, (11) 



where Er = J^— is the recoil energy of the nucleus. The time-dependent function 
T(Er) integrates over all possible kinematic configurations in the scattering process 

r fx 



T{E R ) = 2- i±dv. (12) 

^ •'"mm V 

For a Maxwellian velocity distribution taking into account the motion of the Sun 
and the Earth, Eq. (TT]), one obtains 



T{E R ) = f 



„ , / Vmin + Ve(t)\ „ , ( V min - V e (t) ~ 

Erf - Erf 



(13) 



with ii E R (m x +m A ) 2 
With V mn - ^l 2m 2 mA 



The actually measured response from a nucleus X hit by a WIMP is only a fraction 
of the recoil energy, E = qxE R . It is also called the electron equivalent energy and 
is determined by the quenching factor qx which is different for each detector material. 
Here we assume numerical values of the quenching factors qNa = 0.30 and qi = 0.09 
as they were reported by the DAMA Collaboration || in agreement with previously 
published results |13| . 
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The expected experimental spectrum per energy bin can then be expressed as 

AR f (E+AE)/ qNa dR Na dE R f (E+AE)/ qi dRr dE R 



AE X ' JE/ qNa dEn AE J E / qi dE R y AE 

where r^ a = jt^+m} = 0-153 and rj = M J^^ M] = 0.847 are the respective mass 
fractions of Na and I in the detector material. 

The formulae presented here are in full agreement with those used by DAMA ||. 
The fact that we neglect the WIMP escape velocity from the halo (estimated in Ref. [Q] 
between 580km/s and 625km/s) leads to a negligible numerical difference. The dif- 
ferential event rate in Eq. ([TT] ) is dominated by the scalar contribution due to WIMP 
scattering off iodine. In our calculation we assume that f p ~ /„ |J and use f p as a 
free parameter which can be directly translated into the WIMP cross section on proton 
a p = ^7—^^ I n this approximation one can re- write Eq. ([□]) as 



n(m p +m x ) 2 ■ 

^0.3 



dR 

dE~ = ^ 



R 



-A 2 T(E R )F 2 (q). 



(15) 



Theoretical prediction for the signal is then calculated using Eqs. ( |Tl~D and (|14"D with 
AE = 1 keV. The predicted rates in each channel are evaluated for t = t p and t = 
t p + T/2, and Eq. (|) is used to extract the predicted values Sq\ and S^ k . These are 
compared with the experimental values given in the Table, as described below. 



4 Results and Conclusions 

Our goal is to examine how the sensitivity of the expected event spectrum will depend 
on the velocities of incident WIMPs. It is clear that these enter in a rather complicated 
way, through Eq. (|T3|), but also through the transferred momentum q on which the form 
factors depend. 

Experimental determinations of the Galactic halo velocity distributions are rather 
uncertain. Assuming a spherically symmetric and isotropic halo with a Maxwellian 
dark matter velocity distribution, it was argued in Ref. that t>o should be of order 
of the local Galactic rotation velocity. Several values have been reported for the latter: 
243±20km/s 0, 222±20km/s |TJ, 228±19km/s 0, although a much larger spread 



220 ± 57km/s was also recently quoted |17]]. To be consistent with the DAMA analysis 



where a fixed value of v = 220km/s was assumed, we take v = 220 ± 20km/s. 

We first illustrate in Figure [1] the effect of the dependence of the expected event 
spectrum on the velocities of incident WIMP. We plot the expected differential event rate 
versus the detected energy E for several choices of Vq. We find a significant dependence 
on vq at both smaller and larger values of the detected energy. The effect is caused by 
the dependence of T on the position of the peak of the WIMP velocity distribution t> . 
At fixed m x and a p , T(E) increases with vq for both sodium and iodine. 

The bending appearing in the spectra at lower values of ^| is caused by the con- 
tribution from WIMPs scattering off Na nuclei becoming larger than that off iodine. 
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1 




E (keV) 



Figure 1: The expected differential event rate spectrum versus the detected energy E 
for m x = 60 GeV, a v = 7 x lCT 6 pb and several choices of v . 

The latter decreases more rapidly with energy due to a smaller quenching factor. It is 
worth noting that in the future, with large enough statistics available, the position of 
the bending, which depends on the WIMP's mass and velocity, could provide additional 
information about the signal. 

In order to compare the experimental spectrum measured by DAMA with theoretical 
WIMP predictions, we introduce a function k defined as 




with the experimental errors of both the time-independent and time-dependent signal 
components o" 0i fc and <7 m ,fc, respectively, serving as weights. (They are defined in the 
Table.) Minimization of this function is then used to look for the ranges of WIMP 
masses and cross sections which would best fit the ranges of So,fc and S m ^ responsible 
for DAMA's possible annual modulation signal. 

We stress that the function k cannot serve as a substitute for a dedicated maximum 
likelihood analysis employed by DAMA where detector efficiencies for each crystal and 
other factors were included to evaluate the background and both the signal components 
5o,fc and S m> k- Only a full experimental analysis can lead to selecting a region in the plane 
(m x ,^cr p ) consistent with the signal. Nevertheless, we believe that k is an adequate tool 
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Figure 2: Contours of the function k — 10 (Eq. [T6| ) for different values of the peak of the 
halo WIMP Maxwellian velocity distribution. Denoted by dots is the 2a region selected 
by DAMA assuming vq = 220 km/s. 



to demonstrate the dependence of a signal region on WIMP velocities. This is presented 
in Figure 2 where we plot contours of k = 10 for vq = 220 km/s (solid), which is the 
value used by DAMA, and for vq = 180 km/s (long dash) and Vq = 260 km/s (short 
dash) corresponding to a 2a error in v . We note that the function k is well-focused 
and that the contour k — 10 reproduces the 2a region of DAMA (the dotted curve in 
Figure 2 adopted from Figure 6 in Ref. remarkably well. Any other choice of k 
below 20 or so would give only somewhat more relaxed contours. 

The region of k = 10 slides quickly to the right and becomes much more elongated 
for Vq decreasing from the value Vq = 220 km/s used in the DAMA analysis Q. We find 
a significant relaxation of the upper limit on m x even for values of vq much closer to 
Vq = 220 km/s. For example, for vq = 200 km/s, the region reaches m x w 140 GeV. For 
Vq larger than 220 km/s the region moves to the left and becomes much more confined. 
Its vertical position decreases only very slowly with increasing velocity. 

It is clear from Figure 2 that the dependence of k on vq is very strong and cannot be 
neglected. This effect is particularly important for smaller WIMP velocities as it allows 
for substantially larger values of the WIMP mass to be considered compatible with the 
possible annual modulation signal from DAMA, as already emphasized in Ref. ||18|| . (A 
crude estimate of the effect was also given in Ref. [llj.) We find that varying Vq within 
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Figure 3: Several contours of (a) the time-independent and (b) time-dependent parts 
of k for vq = 220 km/s. Denoted by dots is the contour of 10 of total k. 



the 2a range around 220 km/s (180 < vq < 260 km/s) leads to increasing the upper 
limit for m x by nearly 60% while on the lower side of m x the relaxation is only about 
25%. The range of the £cr p values increases by about 10%. 

Next we examine the relative contributions from the time-independent and time- 
dependent parts to the shape of k contours. This is shown in Figure 3 for the case 
vq = 220 km/s. We find an interesting interplay. The time-independent part of k 
quickly increases around 1.1 x 10 _5 pb (except for smaller m x < 50 GeV) but does not 
put any restriction on large m x . It is the time-dependent part which cuts off the large 
WIMP mass region allowed by the constant component while at the same time favoring 
larger values of £cr p at larger m x . With more data taken over a full period or more 
available, this property may allow one to derive a stronger upper bound on m x , despite 
large uncertainties in halo WIMP velocities. 



In this Letter, we emphasized the importance of astrophysical uncertainties in de- 
termining the shape and overall position of the expected detection spectrum. We did 
this by demonstrating how varying the maximum of the Maxwellian velocity distribu- 
tion of Galactic halo WIMPs will considerably relax the region selected by DAMA as 
corresponding to a possible annual modulation signal. (It is our understanding that the 
DAMA Collaboration is currently in the process of incorporating the effect in their an- 
nual modulation analysis ||20|| .) The effect will also affect WIMP exclusion plots which 
are often drawn for one fixed value of vq, and should be of interest to those involved in 
searches for halo WIMPs. 
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Clearly, there are several other uncertainties entering the halo WIMP detection 
analysis. We reiterate that, while there is much convincing astrophysical evidence for 
the existence of the Galactic halo, parameters describing its shape have only been 
established through indirect methods and their values should be considered as rather 
uncertain. There still remains a considerable spread among possible models for the 
Galactic halo and (halo model dependent) quoted error bars for t> , and the local halo 



density should, in our opinion, be taken with extreme caution. (In Ref. |17| several halo 
models were compared with the standard non-rotating spherical model and their effect 
on the direct rates was found to be minimal.) Nuclear physics uncertainties in the form 



factors are also rather poorly known and have been of much concern |21 . 

We conclude that, due to astrophysical input uncertainties, caution should be taken 
when considering exclusion regions from detection experiments and also in analyzing 
their implications on underlying WIMP models. For example, it has been claimed 



23| that, in supersymmetric models with the lightest neutralino as the stable superpart- 
ner and a WIMP candidate, it is possible to find SUSY configurations consistent with 
the possible annual modulation signal of DAMA but only for large enough tan (3 > 10, 
typically rather small relic abundance, and for restricted ranges of other SUSY parame- 
ters. Including the effect of astrophysical uncertainties will lead to much broader ranges 



of allowed parameters than previously thought ]24 | 
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